Code
# Load libraries
library(terra)
library(sf)
library(tidyverse)
library(tmap)
library(raster)
library(stars)
library(janitor)
library(kableExtra)Leilanie Rubinstein
December 9, 2024
# Import night lights data
night_lights1 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization","data", "VNP46A1/VNP46A1.A2021038.h08v05.001.2021039064328.tif"))
night_lights2 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization", "data", "VNP46A1/VNP46A1.A2021038.h08v06.001.2021039064329.tif"))
night_lights3 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization", "data", "VNP46A1/VNP46A1.A2021047.h08v05.001.2021048091106.tif"))
night_lights4 <- terra::rast(here::here("posts/2024-12-09-houston-blackout-visualization", "data", "VNP46A1/VNP46A1.A2021047.h08v06.001.2021048091105.tif"))# Merge and process night lights data, cropping to the Houston area
houston_extent <- extent(c(-96.5, -94.5, 29, 30.5))
# Check if CRS matches before merging
if (crs(night_lights1) != crs(night_lights2)) {
stop("CRS does not match between night_lights1 and night_lights2")
}
if (crs(night_lights3) != crs(night_lights4)) {
stop("CRS does not match between night_lights3 and night_lights4")
}
night_lights_before <- terra::merge(night_lights1, night_lights2) %>%
crop(houston_extent)
night_lights_after <- terra::merge(night_lights3, night_lights4) %>%
crop(houston_extent)
# Create blackout mask
night_lights_change <- night_lights_after - night_lights_before
night_lights_change[night_lights_change > -200] <- NA
# Vectorize blackout mask
night_lights_change_poly <- as.polygons(night_lights_change) %>%
st_as_sf() %>%
st_make_valid() %>%
st_transform(crs = "EPSG:3083")# Import infrastructure data
roads <- read_sf(here::here("posts/2024-12-09-houston-blackout-visualization",
"data", "gis_osm_roads_free_1.gpkg"),
query = "SELECT * FROM gis_osm_roads_free_1 WHERE fclass='motorway'") %>%
st_transform(crs = "EPSG:3083")
houses <- read_sf(here::here("posts/2024-12-09-houston-blackout-visualization",
"data", "gis_osm_buildings_a_free_1.gpkg"),
query = "SELECT * FROM gis_osm_buildings_a_free_1
WHERE (type IS NULL AND name IS NULL)
OR type in ('residential', 'apartments', 'house', 'static_caravan', 'detached')") %>%
st_transform(crs = "EPSG:3083")# Calculate common min and max for both rasters after log1p transformation
min_val <- min(c(values(log1p(night_lights_before)),
values(log1p(night_lights_after))), na.rm = TRUE)
max_val <- max(c(values(log1p(night_lights_before)),
values(log1p(night_lights_after))), na.rm = TRUE)
par(mfrow = c(1,2))
# Plot before
terra::plot(log1p(night_lights_before),
col = viridis::viridis(50),
main = "Night Lights Intensity \n(February 7, 2021)\n",
zlim = c(min_val, max_val),
plg = list(title = "log1p(nW/cm²/sr)"))
# Plot after
terra::plot(log1p(night_lights_after),
col = viridis::viridis(100),
main = "Night Lights Intensity \n(February 16, 2021)\n",
zlim = c(min_val, max_val),
plg = list(title = "log1p(nW/cm²/sr)"))[1] 157970
houses_count <- houses_filtered %>%
group_by(type) %>%
summarize(count = n()) %>%
ungroup() %>%
st_drop_geometry()
# Check that the sum of buildings adds up to the number of rows in `houses_filtered`
testthat::expect_equal(sum(houses_count$count), nrow(houses_filtered))
houses_count_table <- houses_count %>%
kbl(caption = "Number of homes that experienced power loss (total = 157967)") %>%
kable_classic(html_font = "Cambria")
houses_count_table| type | count |
|---|---|
| apartments | 1136 |
| detached | 353 |
| house | 19760 |
| residential | 1395 |
| static_caravan | 80 |
| NA | 135246 |
Reading layer `ACS_2019_5YR_TRACT_48_TEXAS' from data source
`/Users/leilanie/git/leirubinstein.github.io/posts/2024-12-09-houston-blackout-visualization/data/ACS_2019_5YR_TRACT_48_TEXAS.gdb'
using driver `OpenFileGDB'
Simple feature collection with 5265 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -106.6456 ymin: 25.83716 xmax: -93.50804 ymax: 36.5007
Geodetic CRS: NAD83
Reading layer `X19_INCOME' from data source
`/Users/leilanie/git/leirubinstein.github.io/posts/2024-12-09-houston-blackout-visualization/data/ACS_2019_5YR_TRACT_48_TEXAS.gdb'
using driver `OpenFileGDB'
# Join the median household income from the previous 12 months to the census tract geometries
acs <- left_join(acs_income, acs_texas, join_by(GEOID == GEOID_Data)) %>%
st_as_sf() %>%
st_make_valid() %>%
st_transform(crs = "EPSG:3083")
# Create a hounding box for the Houston area
houston_bbox <- st_bbox(c(xmin = -96.5, xmax = -94.5,
ymin = 29, ymax = 30.5)) %>%
st_as_sfc() %>%
st_set_crs(4326) %>%
st_transform(st_crs(acs))
# Filter to Houston area census tracts
houston_tracts <- acs[houston_bbox, ]
# Identify census tracts that that contained homes that experienced blackouts
houston_tracts_blackouts <- st_filter(houston_tracts, houses_filtered)
# Plot tracts
tm_shape(houston_tracts) +
tm_polygons(col = "white") +
tm_shape(houston_tracts_blackouts) +
tm_polygons(col = "#CBC3E3") +
tm_borders(col = "black") +
tm_compass(type = "8star",
size = 0.8) +
tm_scale_bar() +
tm_layout(main.title = "Census Tracts in Houston that Experienced Blackouts \n(shown in purple)",
scale = 0.8)# Find census tracts that did NOT experience blackouts
houston_non_blackout <- houston_tracts %>%
filter(!GEOID %in% houston_tracts_blackouts$GEOID)
# Label tracts by power status
houston_tracts_blackouts$status <- "Lost Power"
houston_non_blackout$status <- "Did Not Lose Power"
combined_tracts <- rbind(houston_tracts_blackouts, houston_non_blackout)
# Plot combined tracts
ggplot(combined_tracts) +
geom_boxplot(aes(x = B19013e1, fill = status)) + # `B19013e1` = Estimate of median household income in the past 12 months (in 2010 inflation-adjusted dollars)
labs(title = "Distribution of Median Household Income by Power Loss Status",
x = "\nMedian Household Income (in 2011 inflation-adjusted dollars)",
fill = "Power Status") +
scale_fill_manual(values = c("Lost Power" = "#CBC3E3",
"Did Not Lose Power" = "lightgreen")) +
theme_bw() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank())I found that median household income was slightly higher in the census tracts that lost power vs. those that did not lose power. Census tracts on the outskirts of the Houston area generally experienced less power loss, but most of the majority of tracts were still affected. According to Busby et al., all socio-economic groups were affected, but the effects of the storm were harder on low-income families, who live in older, poorly insulated homes and have limited access to resources for repairs and physical health after the storm.
| Data | Citation | Link |
|---|---|---|
| Night Lights | NASA Earth Data (2024). Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC) | Level-1 and Atmospheric Archive & Distribution System Distributed Active Archive Center (LAADS DAAC). |
| Infrastructure Data | Planet OSM (2024). Retrieved from https://download.geofabrik.de/ | Link |
| Census Data | U.S. Census Bureau. (2019). American Community Survey 1-year Public Use Microdata Samples. Retrieved from https://www.census.gov/programs-surveys/acs/news/data-releases.2019.html#list-tab-1133175109 | [Link] (https://www.census.gov/programs-surveys/acs/news/data-releases.2019.html#list-tab-1133175109) |
| Understanding the 2021 Texas Blackout | Busby, J. W., Baker, K., Bazilian, M. D., Gilbert, A. Q., Grubert, E., Rai, V., Rhodes, J. D., Shidore, S., Smith, C. A., & Webber, M. E. (2021). Cascading risks: Understanding the 2021 winter blackout in Texas. Energy Research & Social Science, 77, 102106. https://doi.org/10.1016/j.erss.2021.102106 | Cascading risks: Understanding the 2021 winter blackout in Texas |